Boundary information inflow enhances correlation in flocking 
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The most conspicuous trait of collective animal behaviour is the emergence of highly ordered 
structures. Less obvious to the eye, but perhaps more profound a signature of self-organization, is 
the presence of long-range spatial correlations. Experimental data on starling flocks in 3d show that 
the exponent ruling the decay of the velocity correlation function, C(r) ~ 1/r 7 , is extremely small, 
7 < 1. This result can neither be explained by equilibrium field theory, nor by off-equilibrium 
theories and simulations of active systems. Here, by means of numerical simulations and theoretical 
calculations, we show that a dynamical field applied to the boundary of a set of Heisemberg spins on 
a 3d lattice, gives rise to a vanishing exponent 7, as in starling flocks. The effect of the dynamical 
field is to create an information inflow from border to bulk that triggers long range spin wave modes, 
thus giving rise to an anomalously long-ranged correlation. The biological origin of this phenomenon 
can be either exogenous - information produced by environmental perturbations is transferred from 
boundary to bulk of the flock - or endogenous - the flock keeps itself in a constant state of dynamical 
excitation that is beneficial to correlation and collective response. 
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PACS numbers: 05.65.+b, 87.18-h, 75.10.Hk, 05.50.+q 



Flocking, the collective motion displayed by large 
groups of birds, is one of the most spectacular examples 
of emergent collective behavior in nature, and it has fas- 
cinated inquiring minds since a long time Statistical 
physicists have tackled the problem via minimal models 
of self propelled particles (SPP) [HQ and hydrodynamic 
continuum theories Such studies showed that flock- 

ing can be interpreted as a spontaneous symmetry break- 
ing phenomenon occurring in a "moving ferromagnetic 
spin system" , a sort of non-equilibrium counterpart of 
the well known Heisenberg model 0- The basic ingredi- 
ents of this description - self propulsion, lack of Galileian 
invariance and momentum conservation, local ferromag- 
netic interactions - define an extremely rich universality 
class, able to describe systems as diverse as vertebrate 
herds Q, bacteria colonies @, driven granular matter 
[Toj | , grasshop per swarms 11 1 and active macromolecules 



in living cells [12j, [13 1 . 



Flocking, however, remains a prominent instance of 
collective animal motion for two reasons. First, it in- 
volves large numbers of individuals, hence justifying a 
statistico-mechanical approach to the problem. Second, 
unlike for most 3d systems, for flocks of starlings (Stur- 
nus vulgaris) we have experimental data 14 , , against 
which theories and models can be tested. The statisti- 
cal analysis of individual positions and velocities has re- 
vealed several unexpected physical features that need to 
be explained. In particular, it was found in [1 61 ] that the 
spatial correlations of the velocity fluctuations in starling 
flocks are anomalously long-ranged. Such correlations are 
hard, if not impossible, to reconcile with the current the- 
ories of flocking. 

Consider a flock of N birds with velocities Vj and ve- 
locity fluctuations 5v, = v, — The two point 



correlation function is defined as, 



C(r,L) 
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where CV 



5vi-5vj and rij is the distance between birds 



i and j. In systems of finite size L, due to the global con- 
straint J2i — 0, the function C(r) has a zero, which 
can be used as a finite-size definition of the correlation 
length £, C[r = £) = 0. In starling flocks it was found 
that £ ~ L, namely the correlation function is scale- free 
[l6j |. In fact, long-range correlations are expected in sys- 
tems where a continuous symmetry is spontaneously bro- 
ken (the direction of motion for flocks, the spin direction 
for Hesemberg-like ferromagnets) . However, in starling 
flocks correlations are very long-ranged. 

We can formalize the statement above in the following 
way. Let us write the finite-size correlation function as, 
C(r,L) — £(L)~ 7 g(r/£(L)), where g is a dimensionless 
scaling function, with g(l) = [ijj- Hence, the derivative 
in 1 of the rescaled correlation function, C(r/£), is given 
by C"(r/£ = 1) ~ £(L)-~< ~ £~ 7 . In [HI it was found 
that C"(l) does not show any significant scaling with L 
(or, equivalently, with £) (Fig. implying 7^0 fljj . 
This fact is surprising, as it implies that in two flocks of 
sizes L and 2L, one gets C(2r,2L) ~ C(r,L). Correla- 
tions are basically not decaying at all. 

This result contrasts with the classical Heisenberg 
model on a 3d lattice, where 7=1, implying C(2r, 2L) ~ 
1/2 C (r, L) [7|. The situation does not improve when one 
considers non-equilibrium flocking theories and models. 
While different scalings are expected in the perpendicular 
and parallel directions with respect to the mean velocity, 
the hydrodynamic approach of [H predicts to leading or- 
der 7 = 6/5 in 3d and 7 = 2/5 in 2d [l9j], a result sup- 
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FIG. 1. (color online) (a) Rescaled correlation functions in 3d 
experiments (a) and 2d flocking models (b). Inset: modulus of 
the derivative of the rescaled correlation function in r/£ = 1 
vs. the correlation length £. In both experiment and model, 
correlations are scale free, £ ~ L (not shown here), (a) Exper- 
imental data from different highly ordered flocks with different 
sizes (from 9.1m to 85.7m) and numbers of birds (from 122 
to 4268). (Reprinted from [l6j]) (b) Topological Vicsek model 
fl(]| ; numerical simulations in the highly ordered regime on a 
2d torus (vo = 0.5, N = 256, . . . , 16384, angular noise with 
amplitude -q = 0.15). Data averaged over 5 ■ 10 6 timesteps. 
The inset fit (red line) has slope 7 = 0.4. 



ported by numerical simulations for the 2d topological 
SPP model introduced in ^ (Fig. []}>). 

Therefore, both in 3d and 2d, SPP models and hydro- 
dynamic theories predict a decay of the correlation func- 
tion that is faster than the equilibrium case, whereas in 
starling flocks one finds a decay dramatically slower than 
the equilibrium case. We conclude that the origin of the 
anomalously slow decay of correlation in starling flocks 
is probably not to be found in the self-propelled nature 
of real birds. The discrepancy between models and the- 
ories on one side, and bird flocks on the other side, is 
troublesome. It has been suggested in [16[ that these 
unusually strong correlations may be responsible for the 
very effective response of flocks to external perturbations, 
most notably predators attacks. If this is true, it means 
that the value of 7 may play a relevant evolutionary role. 
Hence, understanding what is going on seems important. 

A first hint about the origin of this phenomenon was 
given in (2lj . There, it was shown that the minimal 
model inferred from the data via a maximum entropy 
criterium, correctly reproduces the slow decay of the cor- 
relation function only if velocities on the flock's boundary 
are kept fixed to their experimental value, while the bulk 
velocities follow the model dynamics. This result sug- 
gests that the slow decay of the correlation function could 
be caused by an information transfer from the boundary 
to the bulk of the flock. However, it is not a priori clear 
which mechanisms is able to enhance correlations in such 
a dramatic way. Here, we hypothesize that this could be 
due to a dynamic information inflow, and test this conjec- 
ture by studying the dynamics of a finite-size Heisenberg 
ferromagnet, under the effect of a dynamical magnetic 
field that affects a part of the boundary. 

We consider the classical Heisenberg model with near- 
est neighbour interactions, defined on a spherical portion 
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FIG. 2. (color online) Time averaged, two points correlation 
of spin fluctuations (r/ — 0.3, a — 2). (a)-(b) Correlation 
function C vs. the rescaled distance r/£ for system radii R 
from 4 to 20 (total number of spins in the legend): (a) weak 
field h = 0.01 and (b) strong field h — 0.5. In the insets: 
correlation length £ vs. R in log-log scale. The dashed red 
lines mark linear growth, (c) Log-log plot of the rescaled 
correlation function slope at r/f = 1 vs. system size R for 
weak (black circles) and strong (red diamonds) fields. The 
dashed blue line marks the decay as 1/R. (d) Magnetization 
m timeseries for the above weak (black line) and strong (red 
line) field cases for R — 15. Inset: Its standard deviation Am 
as a function of field strength. 



(of radius R) of a 3d cubic lattice. The boundary B of the 
sphere (defined as the set of spins with less than 6 neigh- 
bors) is affected by a dynamical external field h* which 
keeps the system out of equilibrium and determines an 
information flow from the boundary to the bulk. The ex- 
ternal field has fixed modulus h, is outward pointing and 
at each time step is applied to only half of the spherical 
boundary. The field dynamics follows a uniform random 
walk on a sphere, designed in such a way that h* reverses 
its direction on average in a time = R a [23j ] . 

Spins s* are unitary vectors with lattice coordinate 
r; = (xi, Tji, Zi), where i = 1, . . . , N and x i} y i7 z, are in- 
tegers such that Ti — ||r,-|| < R. Spins follow the time 
discrete, synchronous dynamics, 



-t+i 



= e 



e 



(2) 



where the sum runs over the lattice nearest neighbors 
of Ti, 6[v] = v/||v|| is a normalization operator and C* 
is a random vector delta correlated in space and time 
and uniformely distributed in the unit spherical sur- 
face. Boundary conditions are determined by the func- 
tion g(rj,h*) = h* if Ti e B and (h* ■ rj) > 0, while 
g(r,-,h*) = otherwise. For zero field, h* = 0, the dy- 
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FIG. 3. (color online) (a) Typical snapshot of spin fluctua- 
tions in the field (red arrow)-magnetization (blue arrow) plane 
for R = 8, h — 0.5, r\ — 0.3, a = 2. The arrows length has 
been rescaled for clarity reasons, (b) Typical snapshot of ve- 
locity fluctuation in a starling flock. (Reprinted from [l6l |) 



namics converges towards the equilibrium distribution of 
an Hcisenberg ferromagnet with a temperature T that is 
a monotonic function of the noise amplitude r\ [22j ]. 

Being interested in the highly ordered phase, we fix 
noise to r] = 0.3. We initially consider a — 2, so that 
the typical field inversion time Th = R 2 is of the order of 
the information propagation time as given by standard 
diffusive dynamics. In order to compare with the results 
of , we define spin fluctuations as u* = s ■ — m*, where 
m' = m*]!* = jj Yli s l is the instantaneous global mag- 
netization and n* its unitary direction. The correlation 
function is defined as in JT]), with Cy = (u- • u*-) f where 
(•)( denotes a time average over a scale r >> r/, [24|. 

The effect of a strong dynamical field on the correlation 
function is striking. In Fig. [2k-b, we report C(r) for 
different 'flock' sizes R, at two values of the field, h = 0.01 
and h = 0.5. In both cases the correlation length £ grows 
linearly with R, as expected in a scale-free system as 
Heisenberg. Moreover, in the weak field case the rescaled 
correlation function C(r/£) behaves as in the equilibrium 
case: the derivative of the correlation function at r/£ = 1, 
vanishes for increasing sizes, C"(l) ~ l/R 1 , with 7=1 
(Fig. black circles). On the contrary, in the strong 
field regime the correlation has a striking resemblance 
with that observed in real bird flocks [l6j . In particular, 
the correlations C(r/£) at different sizes collapse onto a 
single curve. This means that in the strong field regime, 
the derivative C'(l) is constant, implying 7^0 (Fig.[5fc, 
red diamonds), in agreement with experiments (FigJTJi). 

Having enhanced correlation, we must make sure we 
have not destroyed order. Hence, let us check the field 
effects on the global magnetization, which is the equiv- 
alent of flock's velocity. In Figj2ji we show the time se- 
ries of the scalar magnetization m for weak (black) and 
strong field (red). While the strong field standard devi- 
ation Am is increased by about a factor 20, the mean 
magnetization is only slightly reduced, so that the ferro- 
magnet remains in the deeply ordered phase. Note also 
that the variance saturates as the field is increased past 
h = 0.5, so that magnetization fluctuations stay finite 
and relatively small even in the strong field regime. 



We can analytically explain our numerical result by us- 
ing the spin- wave approximation. Let us start by consid- 
ering the Heisenberg model at equilibrium. Each spin can 
be decomposed as Sj = sf n+TTi, where sf and 71^ are the 
longitudinal and perpendicular components with respect 
to magnetization. At low temperature, when the system 
is highly polarized, one has n 2 -C 1, so that, using the 
unitary condition ||s,-|| = 1, we get, sf ~ 1— 7rf/2 (we also 
note that = 7r, at leading order). Under these condi- 
tions the original Hamiltonian, T~L = — 1/2^ « Sj • Sj, 
can be expanded, leading to a Gaussian partition func- 
tion, 



(3) 



where, A, 



• 1 j ■ 



is the discrete Laplacian, and 



the adjacency matrix n,j is 1 for nearest neighbors, and 
otherwise. The connected correlation function Cy can 
be easily computed in terms of the eigenvalues {A a } and 



eigenvectors {w a } (a = 1 • • • N) of Aij [2l|, 
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a>l 



/3A a 



(4) 



where wf = (w a |i), in Dirac's notation. The matrix A 
has a zero eigenvalue related to the original rotational 
symmetry of the Hamiltonian. The first non-zero eigen- 
value is of order l/R 2 on a discrete lattice of size R; it is 
indeed the presence of this soft (or massless) mode that 
gives rise to long-range correlations when R — > 00 [7| . 

Let us first study the effect of a static field, H — > 
H — J2i s i ' hj. The field breaks the rotational symmetry 
and perturbs the diagonal part of the Laplacian matrix, 
An — > Au(h) = J2k Aik + hi ■ n. To first order we have, 
A a ->■ \ a (h) = A a (0) + EiK) 2h i ' n - If th e field is 
homogenous and acts over 0(R 3 ) sites, this correction is 
O(l): eigenvalues are shifted by a nonzero mass and cor- 
relations are no longer scale-free. On the other hand, if 
the field only acts on a number of sites of 0(R 2 ) (as the 
boundary), then we get, X a (h) ~ A a (0) + 0(1/ R), and 
correlations remain long-range. Hence, the first thing we 
learn is that the field must not be applied on all spins 
lest correlation becomes short range. 

From ((3]) we also learn that the correlation function is 
a superposition of normal modes w a (the so-called spin 
waves). Each one of these modes has a specific space 
modulation (on a cubic lattice they are plane waves). The 
lowest non-zero modes correspond to fluctuations with 
length-scale R that reverse the orientation of the spins 
from one side of the system to the other, very similar to 
the fluctuations observed in real flocks (Fig [3]). Hence, 
to lower the value of 7 we must apply a (boundary) field 
that overweights these long-range modes. Since a static 
field leads to 7 = 1 at best, it seems natural to consider 
a time-dependent boundary field, h(t). 
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FIG. 4. (color online) (a) Correlation functions for different 
external field strength (increasing in the blue arrow direction) 
at fixed system size R = 15 and a — 2. In the inset, rescaled 
correlation function slope at r/£ = 1 vs. h in a semi-log scale, 
(b) Correlation functions for different field dynamical expo- 
nent a (increasing in the blue arrow direction) at fixed system 
size R = 15 and field strength h — 0.3. Full symbols refer to 
the "diffusive" exponent a = 2. In the inset, magnetization 
standard deviation as a function of a. In both panels rj = 0.3. 



To treat this case we consider the Langevin equation 
ds, 
It 



j 



hj - /;,s, + 



(5) 



where \ii is a time-dependent Lagrange multiplier en- 
forcing the constraint ||s,|| = 1, and £ 4 is a vectorial 
delta-correlated noise. As in our numerical simulation, 
we choose the subset B where the field is applied to be a 
fraction of the boundary, h»(i) = h for i € B(t), with h 
directed outward. The field has a timescale tk ~ R a ■ If 
temperature is low and the field varies slowly in time, we 
can describe the system in terms of the polarization di- 
rection n(t) and of the instantaneous perpendicular fluc- 
tuations {-Ki(t)} around it. Projecting Eq.([5} along and 
perpendicularly to n, and exploting the equation for the 
constraint we get, 

' /7r; ^ ' -(V^TT.+hH^ 



dt 



M— 

dt 



j 



nf.dn , dn. 



(h • n) ^ 



N 



I^Ch-irO^ + l^ (7) 



where B is the cardinality of B, hi = hfn + hj- and the 
perpendicular component of the noise £ 4 has variance 
4T. From these equations we can gain information on 
how the fluctuations behave. If we assume that n does 
not change significantly on the scale over which fluctua- 
tions decay, then the leading behaviour of the correlation 
function can be estimated using the first line of Eq.^BJ). 
One gets, 

) fce -(A. + A,) t + (5afc ^: (l _ e -2A at) 
a, 6 

+2 f rft'e-^+^X*-*') [m Q (i') ■ hi(t')} } + 6y , (8) 



C i3 {t) 



where b{j is a renormalization term enforcing the con- 
straint, the variables with index a indicate a projection 
on the a normal mode (e.g. h a = J2i w fhi), and the 
magnetization m a (t) is given by, 



m a (i) 



A " f + / dt'e-^-^h^t') . (9) 



The first term on the r.h.s. of Eq.([8]) depends on the ini- 
tial conditions; the second term is the dynamical counter- 
part of the standard Heisenberg correlation and would 
be present even in absence of any external field; the last 
term, on the contrary, is the one mainly affected by the 
presence of a field and by its dynamics. We can see that 
each mode w° gives a contribution to the correlation de- 
caying on a time r a ~ 1/A Q . If the field time-scale, th, is 
much larger than the maximum r a (slowest mode), the 
field is as good as constant, and fluctuations equilibrate 
to their static expression. In this case, correlations are 
of the Heisenberg kind, 7 = 1. 

However, if Th is in the same range as the spin wave 
time scales, all the modes with r a > Th do not equilibrate. 
In particular, we recall that the slowest modes has t q ~ 
1/A a ~ R 2 , and therefore by choosing Th ~ R 2 , we make 
the faster modes relax, but we keep the lowest modes 
excited. The last contribution in Eq.® is therefore non- 
trivial and - if stronger than the second standard term - 
it can modify the large scale behaviour of the correlation. 
For this to occur, we need h 2 ^> f3A a . Given that h a = 

T,ieB w i h i and that ( w i) 2 ~ l l N ( w ° haS n ° rm X )' 
we obtain that the 'strong-field' regime is defined by the 

condition h 2 B 2 /N > f3/R 2 , i.e. TR 3 h 2 > 1. 

This last result tells us that the definition of strong 
vs weak field depends on the size R. For example, for 
h = 0.01 (black circles in Fig|2b) there should be no 
relevant effect of the dynamical field for R < 10 (weak 
field), while some field-induced departure from 7=1 
should be visible for R > 10 (strong field), which is indeed 
what we see in Figl2]:. On the contrary, for h — 0.5 (red 
diamonds in Fij^fc), 7 = 1 is violated at as low a size as 
R ~ 1. In FigfJJi we explicitly show the effect of crossing 
over from weak to strong field. 

We finally checked numerically the effect of changing 
the field timescale, Th (FigHJa). Slow fields (a > 2) pro- 
duce correlation functions very close to the equilibrium 
case, while for a < 2 correlations approach the linear 
decay observed in bird flocks. Therefore, a random dy- 
namical field evolving on a scale faster than R 2 enhances 
correlations in our spin system. However, the effect of 
such fast fields can be detrimental to global order. In the 
inset of FigfJJ} it is shown that the standard deviation 
of the magnetization increases sharply as a —> 1, so that 
fields evolving too fast effectively destroy order in the 
system, whereas dynamics near the timescale, Th ~ R 2 , 
preserves order and increases correlations. 

We note that the optimal timescale, Th ~ R 2 , ex- 
actly characterizes diffusive information propagation in 



5 



the Heisenberg lattice. However, in real flocks birds seem 
to move in the center of mass reference frame in a way 



closer to ballistic than to diffusive, Sr 



1.7 



It 



is therefore possible that the 'right' timescale for enhanc- 
ing the correlation in natural flocks should be somewhat 
faster (i.e. a < 2) than the purely diffusive one. 

Velocity fluctuation correlations have been recently 
measured also in 2d colonies of motile bacteria {Bacil- 
lus subtilis) [27j] ■ This study finds scale- free correla- 
tions, as in the case of bird flocks. However, at vari- 
ance with flocks, the decay exponent found for bacteria 
(7 6 [0.1, 0.2]), is not in plain contradiction with existent 
theories. As we said, the hydrodynamic approach in 2d 
gives 7 = 2/5 4], not far from 0.2, considering exper- 
imental error. Moreover, the linearized hydrodynamics 
theory, which is expected to describe small clusters, and 
thus to be more suited to the data of [27[ , predicts 7 = 



in 2d Q . It therefore seems that in the case of bacteria, 
data can be explained without the need of the border 
perturbation theory developed here. 

We have shown that a dynamical information flow due 
to a fluctuating field and propagating from boundary to 
bulk, gives rise to strong correlations akin to the ones 
observed in real flocks. What is the biological origin (if 
any) of such information inflow? The border of a flock 
is hit by an ever changing flux of environmental stimuli 
and perturbations: attacking falcons, disturbing seagulls, 
wind gusts, sight of significative landmarks, are just few 
examples. These stimuli could give rise to the strong ob- 
served correlation. A few caveats arise about this exoge- 
nous hypothesis. It is unclear whether external stimuli 
can yield a time scale at least comparable to the optimal 
one, Th{R), for every biologically available flock size R. 
To check this point it would be important to study the 
effect of a perturbing field characterized by several differ- 
ent time scales, and see whether the optimal time-scale 
that enhances correlation is naturally selected by the sys- 
tem. Another issue is that there are times at which flocks 
seem not to be subject to evident dynamical perturba- 
tions, and yet display the same anomalous correlations. 

An alternative hypothesis is that the origin of the phe- 
nomenon is endogenous: even in absence of environmen- 
tal perturbations, the flock sets itself constantly into a 
state of dynamical excitation because this behaviour en- 
hances correlation and collective response when true per- 
turbation strikes. Within this scenario, correlation is the 
evolutionary cause of the dynamical excitation, not the 
by-product of it. Moreover, the time scale of the endoge- 
nous excitation would be naturally related to the flock's 
size R, as it is the very birds on the border that spon- 
taneously create the excitation. The problem with this 
endogenous hypothesis is how the flock would do that. 
New models, able to account for spontaneous states of 
dynamical excitation, may help understanding this point. 

This work was supported by grants IIT Seed 
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